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i— I | Abstract 

In the first part of the paper we present a short review of appli- 
r^. ■ cations of digital differential analyzers (DDA) to generation of circles 

0^ \ showing that they can be treated as one-step numerical schemes. In 

^. ' the second part we present and discuss a novel fast algorithm based on 

a two-step numerical scheme (explicit midpoint rule). Although our 
.yv , algorithm is as cheap as the simplest one-step DDA algoritm (and can 

7— ( | be represented in terms of shifts and additions), it generates circles 

£> ■ with maximal accuracy, i.e., it is exact up to round-off errors. 

X: 

Key words and phrases: circle generation, digital differential analyzers, 
exact discretization, explicit midpoint rule 

1 Introduction 

In spite of apparent simplicity fast generation of circles and other curves has 
always been a subject of numerous theoretical and practical studies, find- 
ing applications, among others, in digital plotting, graphical display and 
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numerical machine tool control, compare [El dU EL71 ETJ. Leaving aside the 
fundamental algoritm of Bresenham [3] and techniques based on spline func- 
tions [18], in this paper we focus on digital differential analyzers (DDA), see 
for instance [121 US]- In section [2] we present brief but exhaustive discus- 
sion of digital differential analyzers from a unified point of view, treating 
them as special one-step difference numerical schemes. Then, in section El 
we introduce and discuss a new DDA scheme, cheap and accurate, based on 
a two-step numerical scheme. 

In terms of the natural arc parameter d the equation of the circle of radius 
r is x = rcosfi, y = — rsmfl. The corresponding differential equation is 

dx dy 

M = - y ' M =X - (1) 

Considering various discretizations of ([T|) we can obtain any DDA algorithm. 

2 One-step numerical approximations 



In this section we present a unified approach to digital differential analysers 
representing them by one-step numerical schemes. Therefore they can be 
represented by a general matrix difference equation of the first order: 
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where a, b, c, d depend on the time step h. In order to approximate the system 
(TXJ) these coefficients have to satisfy 
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In other words, A ■■ 
a= l + 0(h 2 ) 



I + hJ + 0(h 2 ), i.e., 
d=l + 0{h 2 ), b = 



-h + 0(h 2 



h + 0(h 2 



The qualitative asymptotic behaviour of the discrete solution x n , y n is deter- 
mined by the characteristic equation 



= det 
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c d — A 



A 2 - (a + d)X + ad-bc 



(4) 



and its roots (eigenvalues of the matrix A) are given by 

Ai = \ (a + d + y/(a + d) 2 - A(ad - be)) = 1 + ih + 0(h 2 ), 
\ 2 =±( a + d- y/(a + d) 2 - 4(ad - be)) = 1 - zft + 0(/i 2 ), 
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2.1 Logarithmic spirals 

The most popular DDA schemes have the form 
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i.e., a — pcosO, c = psin6. Then, 
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Denoting x = r cos <p , y = r Q sin ip , we finally get 
r p n cos(n6 + (p ) , y n = r p n sm(n6 + ip ) . 
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Corollary 2.1. All points generated by ^ lie on the logarithmic spiral de- 
fined by: 
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k:=6- l \np , 



(9) 



where r, ip are polar coordinates on the plane (x,y), compare /!]/. 



From computational point of view the polynomial form of a, c is preferred. 
DDA algorithms can be classified by the order of these polynomials [T2] . 
Below we present several schemes of this kind with the corresponding values 
of p 2 and k (the exact algorithm should have p = 1 and k = 0). 



First order simultaneous DDA algorithm 

a = 1 , c = h, p 2 = l + h 2 , k = ^ - — + • • • ■ (10) 

This is the conventional DDA method |12|. 



Second order simultaneous DDA algorithm 

a=\-\h\ c=h, p 2 = 1 + ^/i 4 , k=^-^ + ..., (11) 
see [HE]. 

Third order simultaneous DDA algorithm 



a = l--h 2 , c =h--h 3 , f ? = l-Lh i + —h«, k = - — + — + ..., 
2 6 P 12 36 24 72 

(12) 
compare [2]. 

Matsushiro's analyzer 

a = l--h 2 , c = h--h 3 , p 2 = l--h A +—h 6 , k = -- + — + . . . . (13) 
2 ' 4 ' ' 4 16 ' 8 48 V ; 

An important point is that this scheme (with h replaced by |/i) was imple- 
mented in terms of shifts and additions (without using multiplications) and 
patented, see [131 E] • 



The best circular interpolator of the third order 

Interestingly enough, the Matsushiro algoritm can be easily improved by 
changing the last term, see [U [T6] : 

a=l--h 2 , c = h--h\ p 2 = l + —h e , fc = — + .... (14) 
2 ' 8 ' ! 64 ' 128 V ; 

We are going to show that this is the best spiral-like circular interpolator of 
the third order in h. 



Lemma 2.2. The most accurate circular interpolator of the form (0|) ; where 
a, c are polynomials of the third order in h, is given by [Lfy 

Proof: We assume 

a = 1 + aih + a 2 h 2 + a 3 h 3 , c = h + c 2 h 2 + c 3 /i 3 (15) 

Then we compute p 2 = a 2 + c 2 : 

p 2 = 1 + 2oi/» + (2a 2 + a 2 + l)/i 2 + 2(a 3 + aia 2 + c 2 )/i 3 + . . . + (a§ + c§)/i 6 . 

Equating to zero five first coefficients we obtain a system of 5 equations for 5 
unknowns, which can be solved easily: 

2oi = => oi = , 

2a 2 + a 2 + l = => 02 = -|, 

«3 + a ia 2 + c 2 = =^> a 3 = -c 2 , (16) 

a\ + 2aia 3 + c 2 + 2c 3 = =^ c 3 = -| - |c| , 

(Z2a 3 + c 2 c 3 = => |c 2 (| - c|) = . 

Then, the coefficient by /i 6 is given by a 3 +c 2 = c 2 + | (| + c 2 ) . The last equation 
of (fT6|) yields two possibilities: either c 2 = or c 2 = |. Corresponding coefficients 
by h e are given by g^ and 1, respectively. Therefore the best accuracy is attained 
in the first case, which leads to (|14p (the second case has another disadvantage: 
coefficients c 2 and a 3 are irrational). 

Interestingly enough, considering /i-expansion of the function k = k(h) we get 
the same final results (although the starting system of equations is more compli- 
cated). The case c 2 = yields k = ^gh 5 + 0(h 6 ) while in the second case, c 2 = |, 
we get k = i/i 5 + 0(/i 6 ). □ 

2.2 Elliptical deformations 

Discrete points generated by d2J) lie on an ellipse if and only if A 2 = Ai and 
| Ai| = 1, which is equivalent to 

\a + d\ < 2 and ad — be = 1 . (17) 

The conditions are known as Barkhausen criteria |19|. 



First order sequential DDA algorithm 
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This is the generator of the so called "magic circle" [131 CEO US] • 

2.3 Elliptical spirals 

Two complex eigenvalues (A2 = Ai) are obtained when 

(a - d) 2 + 46c < , 



(19) 



and, for ad — bc^ 1 the generated circle deforms into an elliptical spiral. The 
following example can be found in [12] . 



Second order sequential DDA algorithm 
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One can easily verify, that (a — d) 2 + 46c = h 2 (h 2 — 4). Therefore f fl9|) is 
satisfied for \h\ < 2. 

Sequentials algorithms are derived from simulataneous algorithms by re- 
placing x n by x n+ i in the equation for y n +i, i-e., simultaneous algorithm (J6j) 
of order M is transformed into 



^ra+l ax n cy n , 

y n+1 = cx n+1 + ay n = acx n + (a - c 2 )y n , 



and then truncated up to the order M, see 
the scheme ( 120]) . 



(21) 



In the case M = 2 we get 



2.4 Exactly circular interpolation 

Digital differential analyzer of the form ([2]) yields an exact circle (up to 
round-off errors) if and only if a — d, c — — b and ad — be = 1. Then, the 
matrix A can be parameterized by a single parameter 9 
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where 4> = (p(h) depends on time step h. In other words, 
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where C 2 (h) + S 2 (h) = 1. Taking 6(h) = h we obtain the exact numerical 
scheme (compare [T71 [20] ) 
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cos h — sin /i 
sin h cos /i 






(24) 



Exact numerical schemes can be constructed for any matrix differential equa- 
tion of the first order, see [5, 6J. 

The Taylor expansion of ( 12"!"]) up to the third order yields the scheme (TT2"|) , 
compare [2]. We recall that this is not the best scheme of the third order, 
see Lemma I 



Implicit midpoint rule 

In the case of linear equations the implicit midpoint rule can be represented 
in an explicit form. Indeed, 



x n +i 



■<- r. 



Vn + Vn+l 



Vn+l ~ Vn 



X n T" 2-n+l 



h 2 ' h 

can be rewritten in a matrix form as follows 
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(26) 



Thus we obtain another algortihm of the form 

4 - h 2 Ah 



where 
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4 + /i 2 



(27) 



see jTTj. Note that a 2 + c 2 = 1 and k = 0, so the implicit midpoint rule 
generates the circle exactly. However, this algorithm has a relatively high cost 
(multiplications in every step). Actually the scheme ( 1271) was a motivation 
for the derivation of the Matsushiro algorithm (TT3|) . see [13]. Indeed, (!T3|) is 
the Taylor truncation of (|27|) up to the third order. 



3 New fast circular interpolator 

The main result of our paper consists in applying a two step method to the 
circular interpolation. The explicit midpoint rule (see [10]), known also as 2- 
step Nystrom method [9] , applied to a general ordinary differential equation 
z = f(t, z) (where zGl™ and / is a given function) reads: 



Z n +2 — Z n + 2hf(t n+ i,Z n+ i) . 

In the case of the circular interpolation z 
Therefore we obtain: 
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compare |16j . This is a two step method and we have to prescribe not only 
xo,yo but also x\,y\. As usual, we can rewrite this scheme in a one-step 
form increasing the dimension of the matrix problem, namely: 
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3.1 First integrals 

In this section we are going to show that ( 129]) preserves circular trajectories 
exactly provided that initial conditions are appropriately chosen. 

Lemma 3.1. If x^ + y^ is a first integral of the system / flPj) . then x n y n+ i — 
x n+iyn is also a first integral. What is more, in this case both first integrals 
are linearly dependent, namely 



X n yn+1 — Xn+iy-n — h{X n + y n ) 



(31) 



Proof: Using (|29|) we compute 

^n+2 + yl+2 = x n + Vn~ ^h(x n y n+ i - x n+1 y n ) + Ah 2 {x 2 n+1 + yl +l ) . 



If x 2 n + yl 



const, then, obviously, x n y n+ i - x n+1 y n 
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Motivated by above results, we define: 
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Lemma 3.2. For any x n ,y n satisfying / TIPj) we have 



(33) 



Proof: A simple straightforward calculation. 



□ 



Theorem 3.3. All points generated by / f^Pj) lie exactly (up to round-off er- 
rors) on the circle x + y = r 2 provided that initial conditions satisfy 



2 _,_ 2 2 



2,2 2 

x 1 + y 1 =r , 



x y 1 - xiy = hr 2 . 



(34) 



Proof: We easily verify that 
Bfi =fx, /i : = 



1 
1 
h 



(35) 



Initial conditions (|46p can be shortly written as £o = f 2 fi- By fj33|) we have 
£n = £o f° r an y n - It means that x 2 n -\- y 2 l = r 2 for any n. □ 



Remark 3.4. The simplest initial conditions satisfying (f^6| ) read 



^o = r , yo = 0, xi = r Vl - h 2 , yi = hr , 
compare FTbJ . 

3.2 Low computational cost 



(36) 



The scheme (129]) has a very low computation cost, to be compared only 
with the simplest one-step scheme (fTOl) . Indeed, both (l29j) and ffTUj) can 
be implemented in terms of two additions and two shifts. However, ( |29|) 
produces an exact circle while the accuracy of ffTUl) is rather limited. 



There is only one more complicated computation, at the start, when 
initial conditions are xi,y% are determined. However, taking into account 
the Taylor expansion 

VT^W = l--h 2 - -h 4 - —h 6 - —h 8 - —h 10 + 0(h 12 ) (37) 
2 8 16 128 256 v ; v ; 

we see that even in that case one can use few shifts to produce x\. Indeed, 
let h = 2~ m and r = 2 N . Then 



Vl-h 2 = 1 - 2~ 2m - 1 - 2" 4m - 3 - 2" 6m " 4 - 2" 8m - 5 - 2~ 8m - 7 + . . . (38) 
and 

_ r)N _ r)N-2m~l _ 9iV-4m-3 _ r)N-6m-4 _ 9^-801-5 _ 9^-8™-? , /on") 

where only integer coefficients are left. 

The case h = ^ (i.e., m = 1) corresponds to the circle approximated by 
regular dodecagon (because then x\ = rcos|, y\ = rsin|, compare (136]) ). 
Assuming iV = 8 we may confine ourselves to the first three terms. 

3.3 Stability 

Characteristic equation det(i? — XI) = reads 

(l-A)(A 2 + 2(l-2/i 2 )A + l) = 0, (40) 

hence eigenvalues of B are given by: 



\ 1 = i, A 2 = 2/i 2 - 1 + 2hVh 2 - 1, A 3 = 2h 2 - 1 - 2hy/h? - 1. (41) 

All the eigenvalues lie on the unit circle: \\k\ = 1 for k = 1, 2, 3. Correspond- 
ing eigenvectors are denoted by fk, k = 1,2,3, namely 

Bfi = h , Bf 2 = A 2 / 2 , Bf 3 = A 3 / 3 , (42) 

(/1 is given by ( !35l) and the form of other eigenvectors is not important). 
The initial condition £ can always be represented as a linear combination 
of eigenvectors: £0 = 0/1 + 6/2 + c/3, where we assume a ~ r 2 (compare 
Theorem 13. 3p and h, c are small perturbations. Using (142]) and (133]) we obtain 

Zn = fi + b\%h + c\ n 3 f 3 . (43) 

Taking into account | A^l = | A3 1 = 1, we see that the initial perturbations 
remain unchanged during the evolution. 

10 



3.4 Exact circular interpolator 

The scheme (1291) produces exact circular trajectory but the period is modified. 
For instance, in the case of regular dodecagon (h = |) the total change of $ 
after 12 steps is A# = 12 h = 6 while the exact value is, of course, 2n. In 
general, the period of the circular interpolator ( 1291) is given by 

2tt/i 

T = — . 44 

arcsin h 



In order to preserve exactly not only the trajectory but also the period of 
the circular motion, we can use a nonstandard modification of ( 1291) . following 
the approach of |5j |6] : 

X n +2 = x n — 2oy n j r \ , 

(45) 

Un+2 = Vn + 25x n+1 , 

where 6 = 6(h) is a given function of h. 

Theorem 3.5. For any 5 = 5(h) all points generated by JHJ5\ ) lie exactly (up 
to round- off' errors) on the circle x + y = r 2 provided that initial conditions 
satisfy 

2,22 2,22 r2 (ac\ 

x o + Vo = r . x i + Vi = r , x y! - x t y = br . (46) 

Proof: It is enough to repeat the proof of Theorem 13.31 (with h replaced by 5 in 
appropriate places). □ 



The period of the circular interpolator ( 1451) is 

2tt/i 

T = — vh\- 47 

arcsin oyti) 



Theorem 3.6. If 5 = sinh, then scheme 
x n +2 = x n - 2y n+1 sin h , 

Vn+2 = Vn + 2x n+ i sin h , (48) 

x yi - Xxy = r 2 sin h , 

generates the circle exactly (up to round-off errors) preserving the period. 

11 



Proof: The scheme (|45p will preserve the period exactly if x n = r cos hn and 
y n = r sin hn. Taking into account well known trigonometric identities 

cos(hn + 2h) — cos(hn) = —2 sin h sin( hn + h) , 

(49) 
sin(/m + 2/t) — sin(/m) = 2 sin /i cos(hn + /i) , 

we conclude that scheme (|45p is exact for 5{h) = svuh. □ 

The exact interpolator (I48p uses two multiplication by sin h at every step, 
hence it is relatively expensive (but two times cheaper than its one-step 
conterpart (JM]))- 

The scheme (145]) can serve as a starting point for deriving cheaper algo- 
rithms by taking /i-polynomials instead of sin h. All these schemes preserve 
exactly the circular trajectory. The best approximation of the period are 
obtained for Taylor truncations. In particular, taking 

S = h - \h* (50) 

6 

we get the best scheme fj4"5]) of order 3. The multiplication by | can be 
replaced by shifts because 



DC 



- = , l ^ = - y 4- fc = y 2- 3 - 2fc , (si) 

and we can approximate 5 by 5n given by 

N 

d N = h-h 3 y^2- 3 ~ 2k , (52) 

k=0 

where iV is chosen to assure the required accuracy. 

4 Conclusions 

We discussed DDA circular interpolators based on one-step numerical schemes. 
We have shown that the best cheap interpolator of the third order is given by 
( I14p . see Lemma \2. 21 This algorithm can be expressed in terms of shifts and 
additions only, and it has quite good accuracy. Exact interpolators based on 
one-step methods are more expensive because multiplications in every step 
are required. 

12 



In section [3] we presented a family of circular interpolators, based on a 
two-step method, namely the explicit midpoint rule. All of them interpolate 
the circle exactly (up to round-off errors). The simplest interpolator, given 
by (129]) . is very cheap and is very good for graphical applications (when the 
period of the circular motion is irrelevant). The scheme ()48p preserves exactly 
also the period but is more expensive. Taking an appropriate polynomial 
5(h) we can control the accuracy and computational cost of the scheme (PR)]) , 
obtaining interpolators suitable for different purposes. 
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